The R Project for Statistical Computing

• https://www.r-project.org
• Software estadístico libre, abierto y gratuito
• Miles de librerías y proyectos complementan el software base en distintas plataformas:
  • CRAN
  • R Forge
  • https://www.bioconductor.org/
• Una publicación:
  • R News (2001 - 2008)
  • R Journal (2009 - presente)
• Una gran comunidad de usuarios escribe en blogs y foros:
  • https://www.r-bloggers.com/
  • https://rweekly.org/
  • https://stackoverflow.com/questions/tagged/r
• Por qué R?
  • Porque podéis escribir las instrucciones en R scripts, y ejecutarlas todas de una vez, sin tener que memorizar y repetir cada uno de los pasos.
  • Porque en Internet podréis encontrar solución a (casi) todos los problemas que os encontréis
  • Porque podréis personalizar/modificar cada instrucción con las opciones deseadas.
  • Porque podéis trabajar con más de una base de datos abierta simultáneamente y no tienen que estar guardadas en un formato específico para poder abrirlas.

R como calculadora

• R es esencialmente una consola en la que el cursor se sitúa tras el símbolo >. Ahí se escriben las instrucciones. Se ejecutan con Enter.
• Los resultados suelen aparecer debajo. En el caso de gráficas, depende del entorno en que se trabaje (R, RStudio, R Commander, etc.)
• Para empezar, R permite realizar las operaciones de cualquier calculadora científico.
• Incluye además instrucciones específicas para operaciones estadísticas.

Necesitamos un editor

• Más práctico: ejecutar múltiples operaciones simultáneamente.
• Para ello: se escriben todas en un fichero de texto simple, pero con extensión .r en lugar de .txt.
• Es suficiente un editor de texto simple. En Windows: Bloc de notas (en OS X: TextEdit). Los editores con formatos complejos como Microsoft Word no valen.
• Otros editores facilitan la introducción y ejecución de las instrucciones.
• Ej.: Notepad++, Sublime Text, Visual Studio Code, Geany, Nano, Emacs, VIM, …
• Enfocado al uso de R y el más usado para ello: RStudio (nombre del editor y de la empresa que lo produce y comercializa).
  • RStudio es una interfaz que tiene el objetivo de facilitar y agilizar el uso de R y otros formatos que interactuan con R (documentos pdf, archivos word, páginas y aplicaciones web, etc., todos ellos con código R ejecutado).
  • En particular ofrece menús para especificar el directorio de trabajo, instalar, actualizar y/o cargar librerías, y abrir y visualizar bases de datos, entre otros. Asímismo dispone de atajos de teclado para introducir los comandos más habituales de manera más rápida.
  • RStudio no es software libre, aunque su edición básica es de código abierto y gratuita.
  • https://www.rstudio.com/

• Permite tener un fichero con extensión .r en el que se puede escribir el código (con autocompletado) y ejecutarlo
• Permite añadir instrucciones mediante atajos de teclado

Algunas nociones sobre RStudio

• Las instrucciones en el R-script en RStudio se ejecutan con Control+Enter

Qué es lo primero que nos encontramos en R? Objetos, funciones, comentarios, (directorios) y librerías

• Objetos elementales:
  • Numéricos, como en una calculadora (vimos algunos hace unas diapositivas)
  • Cadenas de texto, siempre limitadas por comillas simples o dobles, desde un carácter, por ejemplo "a", o 'z', hasta un texto largo como "Lorem ipsum dolor sit amet..." o 'Lorem ipsum dolor sit amet...'
• Funciones: ejecutan las instrucciones que las definen. Genéricamente se ejecutan escribiendo el nombre de una función f seguido de los objetos o valores que necesita entre paréntesis.
  • Cuando no es necesario ningún valor, estos paréntesis se escriben sin incluír nada, así: f()
  • En la mayoría de casos, necesitan evaluar objetos. La evaluación de un objeto x con una función f se codifica f(x).
• Un comentario es un fragmento de código, texto… no ejecutable. Se señalan situando a la izquierda un símbolo #
• Directorio de trabajo: el entorno de R tiene como referencia un directorio concreto. Esto nos permite no tener que introducir la dirección de la carpeta cuando queremos acceder a algún elemento de esta (por ejemplo, ficheros con bases de datos).
  • Para saber cuál es el directorio de trabajo actual utilizamos la función getwd;
  • Para cambiar el directorio se utiliza la función setwd evaluando el directorio de trabajo en que queremos situarnos ahora;
  • Se puede saber qué archivos hay en el directorio donde R está trabajando con la función list.files;

getwd() # el directorio donde R está trabajando
setwd("P:/projects/Curso-R") # (en Windows) se tienen que sustituir los \ por /
list.files()

O bien, mediante RStudio:

`Session > Set Working Directory > Choose Directory`
`Session > Set Working Directory > To Source File Location`

• Librerías:
  • Para instalar una librería llamada tidyverse evaluamos con la función install.packages la cadena consistente en el nombre de la librería, "tidyverse";
  • Para cargar una librería llamada readxl, evaluamos con la función library (o bien require, una función distinta), bien la cadena con el nombre de la librería "readxl" o bien (el nombre de) la librería misma, readxl.

install.packages("tidyverse") # instala la librería tidyverse, que a su vez instala las librerías readxl, dplyr, tidyr, haven, y otras (https://www.tidyverse.org, https://github.com/tidyverse/tidyverse)
library(readxl)
# library("readxl")
# require(readxl)
# require("readxl")

O bien, mediante RStudio:

Panel `Packages > Install`

Siguientes pasos: formas de obtener ayuda, argumentos de una función y abrir un fichero con formato de Excel

• R dispone de una ayuda muy completa sobre todas las funciones, procedimientos y elementos que configuran el lenguaje
• Se puede acceder a ella con la función help() o mediante ?. Los siguientes ejemplos nos muestran los documentos de ayuda de las respectivas funciones.

help(install.packages)
?library

• Para evaluar un objeto mediante una función, implícitamente estamos asignando el objeto a un argumento de la función. Por ejemplo, al ejecutar library(readxl) asignamos (el nombre de) la librería readxl al primer argumento de la función library, que es package. Entonces, otra forma de ejecutar library(readxl) es ejecutar library(package = readxl). En este caso, es exactamente lo mismo, pero cuando queremos usar otro argumento de una función distinto del primero en general tenemos que nombrarlo.
• Además de las funciones, entre otros objetos de los que podemos ver documentos de ayuda están las librerías.

?help
help(package = "readxl") # Informémonos acerca de la librería readxl
?read_excel # ayuda de la función read_excel

Para abrir ficheros de excel disponemos de diversas funciones en las librerías readxl y openxlsx entre muchas otras. Aunque la segunda librería tiene más funciones para manipular y escribir ficheros .xlsx, la primera dispone de funciones que nos permiten leer también ficheros .xls

read_excel("data_curs_stat/EP1.xls")

## # A tibble: 4,753 × 12
##    q0002_hhid sex   age   mar_1 edu1      phys_hea1                    hea1 dep1  score_lon1 score_sup1 income1 income_inf1
##         <dbl> <chr> <chr> <chr> <chr>     <chr>                       <dbl> <chr>      <dbl>      <dbl>   <dbl> <chr>      
##  1          1 fem   35-49 Yes   Tertiary  None                         70.5 No             3         12      15 Normal     
##  2          2 masc  35-49 Yes   Tertiary  None                         78.9 No             3         12      17 Normal     
##  3          3 masc  65-79 Yes   Primary   None                         66.6 No             3         12      15 Normal     
##  4          4 masc  50-64 Yes   Primary   1 physical health problem    77.4 No             3         11      18 Normal     
##  5          5 fem   50-64 No    Tertiary  None                         52.5 No             3         12      29 Normal     
##  6          6 fem   50-64 No    Tertiary  None                         53.7 No             3         13      25 Normal     
##  7          7 masc  50-64 Yes   Secondary None                         59.1 No             6         10      12 Good       
##  8          8 masc  50-64 Yes   Primary   None                         76.4 No             3         12      14 Normal     
##  9          9 fem   50-64 No    Secondary None                         73.7 No             3         11      26 Good       
## 10         10 fem   80+   No    Primary   2+ physical health problems  28.0 No             6         13      10 Normal     
## # … with 4,743 more rows

Abrir bases de datos y ejecutar funciones: R no es como Stata.

La función read_excel abre la base de datos, pero a diferencia de instrucciones como use o import excel en Stata, no la copia a memoria. Así, si en Stata podríamos ejecutar describe o summarize, si en R ejecutamos una función como summary el resultado es la propia función.

summary

## function (object, ...) 
## UseMethod("summary")
## <bytecode: 0x55d5dab07f10>
## <environment: namespace:base>

Tenemos que evaluar la base de datos que abrimos con read_excel con la función summary:

summary(read_excel("data_curs_stat/EP1.xls"))

##    q0002_hhid        sex                age               mar_1               edu1            phys_hea1              hea1           dep1             score_lon1     score_sup1   
##  Min.   :    1   Length:4753        Length:4753        Length:4753        Length:4753        Length:4753        Min.   : 0.00   Length:4753        Min.   :3.00   Min.   : 3.00  
##  1st Qu.: 1327   Class :character   Class :character   Class :character   Class :character   Class :character   1st Qu.:42.83   Class :character   1st Qu.:3.00   1st Qu.:10.00  
##  Median : 2615   Mode  :character   Mode  :character   Mode  :character   Mode  :character   Mode  :character   Median :55.64   Mode  :character   Median :3.00   Median :12.00  
##  Mean   : 3145                                                                                                  Mean   :53.97                      Mean   :3.74   Mean   :11.53  
##  3rd Qu.: 3807                                                                                                  3rd Qu.:66.06                      3rd Qu.:4.00   3rd Qu.:13.00  
##  Max.   :92755                                                                                                  Max.   :92.82                      Max.   :9.00   Max.   :14.00  
##                                                                                                                 NA's   :187                        NA's   :285    NA's   :414    
##     income1       income_inf1       
##  Min.   : 1.000   Length:4753       
##  1st Qu.: 2.000   Class :character  
##  Median : 5.000   Mode  :character  
##  Mean   : 8.989                     
##  3rd Qu.:15.000                     
##  Max.   :35.000                     
##  NA's   :1067

Abrir la base de datos mediante la función read_excel cada vez que queremos operar con ella no es práctico. En lugar de esto, guardamos la base de datos en una variable.

Objetos complejos. Bases de datos

R es un lenguaje de programación que nos permite acceder a objetos con distintas estructuras y manipularlos. Los tipos de objetos más importantes son:

• El vector, una colección de datos elementales del mismo tipo básico ordenados.
  • Los tipos básicos más comunes son:
    • Numérico
    • Cadena de texto o caracteres
    • Lógico: TRUE, FALSE y NA (missing)
  • Habitualmente, los vectores se crean y se unen con la función c
• La lista, una colección de objetos de cualquier tipo. Habitualmente se crean mediante la función list
• La función
• El símbolo, el nombre de un objeto
• El objeto nulo o vacío, NULL

La mayoría de los objetos pueden tener atributos, como pueden ser las dimensiones o las etiquetas. De esta manera se pueden crear objetos con estructuras más complejas:

• Matrices (resp. arrays), son vectores ligados a un atributo dimincluyendo las 2 (resp. -) dimensiones del objeto.
• Factores, vectores con un atributo levels que corresponde a las posibles categorías de una variable categórica y un atributo class="factor".
• Data frames, el tipo básico de base de datos en R. Es una lista de vectores (las columnas de la base de datos) -que pueden ser matrices o factores- que tienen todos la misma longitud, con un atributo names que designa los nombres de las columnas o variables de la base de datos.

Disponemos de un conjunto de operadores en R que nos permiten acceder y manipular los objetos. Entre ellos:

• <-, para asignar un nombre a un objeto. Por ejemplo, x <- c(3,log(4), 5.2), y <- list(3, 'a', list(TRUE))
• [, para acceder a un elemento de un vector o a una sublista. Por ejemplo x[3], y[2]
• [[, para acceder a un elemento en una lista. Por ejemplo y[[2]]

Podemos informarnos acerca de la clase, estructura y atributos de un objeto por medio de las funciones class, str y attributes respectivamente.

• Los missings en R se representan con el símbolo NA, existente para los distintos tipos básicos, o también NaN en el caso de objetos numéricos.

Primeros pasos con una base de datos.

Asignamos la base de datos a una variable que llamamos dataw1:

dataw1 <- read_excel("data_curs_stat/EP1.xls")
#summary(dataw1)
#str(dataw1)
class(dataw1)

## [1] "tbl_df"     "tbl"        "data.frame"

head(dataw1) # primeras (6) filas

## # A tibble: 6 × 12
##   q0002_hhid sex   age   mar_1 edu1     phys_hea1                  hea1 dep1  score_lon1 score_sup1 income1 income_inf1
##        <dbl> <chr> <chr> <chr> <chr>    <chr>                     <dbl> <chr>      <dbl>      <dbl>   <dbl> <chr>      
## 1          1 fem   35-49 Yes   Tertiary None                       70.5 No             3         12      15 Normal     
## 2          2 masc  35-49 Yes   Tertiary None                       78.9 No             3         12      17 Normal     
## 3          3 masc  65-79 Yes   Primary  None                       66.6 No             3         12      15 Normal     
## 4          4 masc  50-64 Yes   Primary  1 physical health problem  77.4 No             3         11      18 Normal     
## 5          5 fem   50-64 No    Tertiary None                       52.5 No             3         12      29 Normal     
## 6          6 fem   50-64 No    Tertiary None                       53.7 No             3         13      25 Normal

#head(dataw1, 3) # primeras 3 filas
#tail(dataw1) # últimas (6) filas

Ejercicio:

Usando las funciones setwd y read_excel convenientemente (hay varias posibilidades), guardar también las bases de datos EP2.xls y EP3.xls en dos variables que llamaremos dataw2 y dataw3. Podéis ver dataw1, dataw2, dataw3 y sus dimensiones en el panel Environment de RStudio?

Nota: Para abrir ficheros *.csv disponemos de la función read.csv en R, de la función read_csv en la librería readr y de la función fread en la librería data.table entre otras. La librería haven dispone de diversas funciones que permiten abrir y escribir ficheros en formatos de Stata, SPSS y SAS (pueden consultarse en su ayuda, package(help = "haven")). Por otra parte la librería readspss (https://github.com/JanMarvin/readspss) tiene funciones que permiten abrir bases de datos en formato de SPSS encriptadas con contraseña.

Acceso a las variables de una base de datos: los operadores $ y [[

# dataw1$phys_hea1
# dataw1[["phys_hea1"]]
class(dataw1$phys_hea1) # chr quiere decir que está guardada como una cadena de carácteres

## [1] "character"

#summary(dataw1$phys_hea1)
str(dataw1$phys_hea1)

##  chr [1:4753] "None" "None" "None" "1 physical health problem" "None" "None" "None" "None" "None" "2+ physical health problems" "None" "None" "2+ physical health problems" ...

#attributes(dataw1$phys_hea1)
table(dataw1$phys_hea1)

## 
##   1 physical health problem 2+ physical health problems                        None 
##                        1538                         887                        1723

#?table
table(dataw1$phys_hea1, useNA = "ifany")

## 
##   1 physical health problem 2+ physical health problems                        None                        <NA> 
##                        1538                         887                        1723                         605

prop.table(table(dataw1$phys_hea1))

## 
##   1 physical health problem 2+ physical health problems                        None 
##                   0.3707811                   0.2138380                   0.4153809

Guardando los resultados

Hemos dicho que R es como una calculadora, y que si no asignamos los objetos a variables, se muestran en consola pero no quedan guardados en ningún sitio.

• Los objetos de R (como una base de datos o dataframe) se pueden guardar con la función save en ficheros con la extensión .rda o .rdata (aunque a veces también se escribe la r e incluso la d en mayúscula, por ejemplo .RData)
• Los resultados, si son matrices, data frames o tablas, además, también se pueden guardar en ficheros .csv o con formatos de excel (a menudo, las funciones mencionadas de la forma read_ o read. tienen correspondientes write_ o write.).
• Los resultados y más generalmente todo aquello que aparece por consola puede ser guardados en ficheros de texto con la función sink.
• A no ser que se especifique la dirección del fichero, serán guardados en el directorio donde R está trabajando (getwd())

#?sink
sink("Summary_Estudi_poblacional_w1.txt", split = TRUE)
table(dataw1$phys_hea1)

## 
##   1 physical health problem 2+ physical health problems                        None 
##                        1538                         887                        1723

summary(dataw1$hea1)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##    0.00   42.83   55.64   53.97   66.06   92.82     187

sink()

Ejercicio:

Qué aparece en el fichero Summary_Estudi_poblacional_w1.txt?

Transformación de variables (I)

• La función as.factor permite transformar a una variable categórica

class(as.factor(dataw1$phys_hea1))

## [1] "factor"

table(as.factor(dataw1$phys_hea1), useNA = "ifany")

## 
##   1 physical health problem 2+ physical health problems                        None                        <NA> 
##                        1538                         887                        1723                         605

• En formatos de fichero distintos de .xls o .csv, como pueden ser los formatos de Stata o SPSS, las categorías de las variables categóricas pueden estar grabadas como etiquetas, mientras los valores son numéricos. En esos casos, para recuperar las etiquetas como categorías es más conveniente usar la función as_factor de la librería haven.

Guardamos la transformación en la base de datos. Para ello usaremos la función mutate de la librería dplyr, que permite realizar varias transformaciones separadas por comas

library(dplyr)
#?mutate

dataw1 <- mutate(dataw1, phys_hea1 = as.factor(phys_hea1))

# es lo mismo que:

dataw1 <- dataw1 %>% mutate(phys_hea1 = as.factor(phys_hea1))

# podemos escrbir el argumento a la derecha de %>% en las líneas inferiores

dataw1 <- dataw1 %>% #partimos de la base de datos dataw1 y entonces
  mutate(phys_hea1 = as.factor(phys_hea1)) # transformamos la variable y la guardamos con el mismo nombre

# finalmente la base de datos queda guardada con el mismo nombre mediante la asignación inicial (dataw1 <- ...)

• El primer argumento o valor de la función mutate es la base de datos o data frame.
• Los siguientes argumentos son las transformaciones variable_output = transformación separadas por comas.
• Usando %>%, la función situada en el lado derecho (en este caso mutate) ya reconoce que su primer argumento es el valor situado en el lado izquierdo de %>% (en este caso dataw1)

Pipe

Transformación de variables (II)

Como la variable q0002_hhid es un id podría interesarnos que fuera de clase cadena

#class(dataw1$q0002_hhid)
#?as.character

Varias transformaciones seguidas las podemos evaluar: - Aplicando varias veces mutate:

dataw1 <- dataw1 %>% #partimos de la base de datos dataw1 y entonces
  mutate(phys_hea1 = as.factor(phys_hea1))  %>% # transformamos la variable phys_hea1 y la guardamos con el mismo nombre, y entonces
  mutate(dep1 = as.factor(dep1)) %>% # transformamos la variable dep1 y la guardamos con el mismo nombre
  mutate(q0002_hhid = as.character(q0002_hhid)) # transformamos la variable q0002_hhid y la guardamos con el mismo nombre

• Pero en lugar de eso, habitualmente se separan las instrucciones por comas en un sólo mutate:

dataw1 <- dataw1 %>% #partimos de la base de datos dataw1 y entonces
  mutate(phys_hea1 = as.factor(phys_hea1), dep1 = as.factor(dep1), q0002_hhid = as.character(q0002_hhid)) # transformamos las variables phys_hea1,  dep1 y q0002_hhid y las guardamos con los mismos nombres

Ejercicio:

Buscar todas las variables categóricas en las 3 bases de datos y transformar estas bases de datos de manera que esas variables sean factores. Mostrar las frecuencias de las categorías de algunas de esas variables. Transformar la variable q0002_hhid a cadena de carácteres. Consultar también la ayuda de as.numeric, as.integer y as.factor

Otras operaciones con funciones de la librería dplyr y funciones de R para el cálculo de estadísticos

• Ejemplo 1

dataw1 %>% #partimos de la base de datos dataw1 y entonces
  select(q0002_hhid, hea1, dep1) %>% # mantenemos sólo las columnas q0002_hhid, hea1 y dep1, y entonces
  filter(dep1 == "Yes") %>% # nos quedamos con las filas de quienes padecen depresión, y entonces
  arrange(desc(hea1)) %>% # ordenamos las filas por de mayor a menor valor de estado de salud, y entonces
  head(3) # nos quedamos con los 3 pacientes con depresión con mayor valor de Health state

## # A tibble: 3 × 3
##   q0002_hhid  hea1 dep1 
##   <chr>      <dbl> <fct>
## 1 4941        82.1 Yes  
## 2 4339        80.2 Yes  
## 3 1299        79.8 Yes

• Ejemplo 2 (uso de funciones de R como mean, median o quantile; se pueden usar otras como var, sd, min, max, IQR, etc.)

dataw1 %>% #partimos de la base de datos dataw1 y entonces
  select(q0002_hhid, hea1, dep1) %>% # mantenemos sólo las columnas number_id, hea1 y dep1, y entonces
  group_by(dep1) %>% # agrupamos por las categorías de depresión, y entonces
  summarise(mean_hea1 = mean(hea1, na.rm = TRUE), median_hea1 = median(hea1, na.rm = TRUE), tercil_hea1 = quantile(hea1, probs = 1/3, na.rm = TRUE), n = n()) # calculamos la media, la mediana y el primer tercil de hea1 para cada categoría de depresión y el número de observaciones por cada categoría

## # A tibble: 3 × 5
##   dep1  mean_hea1 median_hea1 tercil_hea1     n
##   <fct>     <dbl>       <dbl>       <dbl> <int>
## 1 No         55.9        57.4        50.2  4073
## 2 Yes        38.4        37.9        31.1   510
## 3 <NA>      NaN          NA          NA     170

Nota: usamos na.rm = TRUE dentro de mean, de median y de quantile para que calcule la media de aquellos valores que no son missing. En caso contrario, cuando hay missings el resultado es NA.

Nota: Todo lo anterior se puede realizar también con funciones de R sin necesidad de acudir a la librería dplyr, pero las alternativas pueden ser más complejas.

Cálculo de prevalencias:

Prevalencia de depresión:

dataw1 %>% 
  count(dep1) %>% # Contamos los individuos en cada categoría de depresión y entonces
  mutate(prop = 100*n/sum(n)) # calculamos el %

## # A tibble: 3 × 3
##   dep1      n  prop
##   <fct> <int> <dbl>
## 1 No     4073 85.7 
## 2 Yes     510 10.7 
## 3 <NA>    170  3.58

Prevalencia de depresión por grupos de edad y sin tener en cuenta los missings:

# ?is.na
dataw1 %>% 
  group_by(age) %>% 
  count(dep1) %>% # Contamos los individuos en cada categoría de depresión por grupo de edad y entonces
  filter(!is.na(dep1)) %>% # eliminamos los missings si nos interesa contar el porcentaje sobre el total de respuestas válidas
  mutate(prop = 100*n/sum(n)) # Calculamos el %

## # A tibble: 10 × 4
## # Groups:   age [5]
##    age   dep1      n  prop
##    <fct> <fct> <int> <dbl>
##  1 18-34 No      382 93.9 
##  2 18-34 Yes      25  6.14
##  3 35-49 No      500 90.7 
##  4 35-49 Yes      51  9.26
##  5 50-64 No     1558 88.5 
##  6 50-64 Yes     202 11.5 
##  7 65-79 No     1298 87.3 
##  8 65-79 Yes     188 12.7 
##  9 80+   No      335 88.4 
## 10 80+   Yes      44 11.6

Variables y dimensiones de las bbdd

names(dataw1)

##  [1] "q0002_hhid"  "sex"         "age"         "mar_1"       "edu1"        "phys_hea1"   "hea1"        "dep1"        "score_lon1"  "score_sup1"  "income1"     "income_inf1"

dim(dataw1)

## [1] 4753   12

ncol(dataw1)

## [1] 12

nrow(dataw1)

## [1] 4753

names(dataw2)

## [1] "q0002_hhid" "dep2"       "score_lon2" "score_sup2" "income2"    "phys_hea2"  "hea2"

dim(dataw2)

## [1] 2400    7

names(dataw3)

## [1] "q0002_hhid" "dep3"       "score_lon3" "score_sup3" "income3"    "phys_hea3"  "health3"

dim(dataw3)

## [1] 1512    7

Fusión (merge) de bbdd

Cuántos missings tiene la variable q0002_hhid?

table(is.na(dataw1$q0002_hhid)); table(is.na(dataw2$q0002_hhid)); table(is.na(dataw3$q0002_hhid))

## 
## FALSE 
##  4753

## 
## FALSE 
##  2400

## 
## FALSE 
##  1512

Ya tenemos las bbdd abiertas y guardadas. Sólo tenemos que unir por la(s) variable(s) identificadora(s), en este caso q0002_hhid.

data <- dataw1 %>% # partimos de la base de datos dataw1 y entonces
  full_join(dataw2, by = "q0002_hhid") %>% # unimos horizontalmente con todas las observaciones de dataw2 con q0002_hhid iguales a los de dataw1 y añadimos las nuevas, y entonces
  full_join(dataw3, by = c("q0002_hhid"))# unimos horizontalmente con todas las observaciones de dataw3 con q0002_hhid iguales a los que ya había y añadimos las nuevas
names(data)

##  [1] "q0002_hhid"  "sex"         "age"         "mar_1"       "edu1"        "phys_hea1"   "hea1"        "dep1"        "score_lon1"  "score_sup1"  "income1"     "income_inf1"
## [13] "dep2"        "score_lon2"  "score_sup2"  "income2"     "phys_hea2"   "hea2"        "dep3"        "score_lon3"  "score_sup3"  "income3"     "phys_hea3"   "health3"

dim(data)

## [1] 4753   24

?full_join

Ejercicio:

Con full_join creamos una base de datos resultado de fusionar las 3 iniciales e incluír todas las observaciones de cada una de ellas. Consultando en la ayuda, este ejercicio consiste en fusionar las 3 bases de datos, pero incluyendo sólo aquellas observaciones de comunes a las 3. Cuántas observaciones tiene?

Combinar bases de datos

Abrimos las 3 olas de un ensayo clínico en fichero separados.

ac1 <- read_excel("data_curs_stat/AC1.xls")
ac2 <- read_excel("data_curs_stat/AC2.xls")
ac3 <- read_excel("data_curs_stat/AC3.xls")

dim(ac1); names(ac1); dim(ac2); names(ac2); dim(ac3); names(ac3)

## [1] 400   7

## [1] "q0002_hhid" "sex"        "age"        "mar_1"      "edu1"       "hea1"       "grups"

## [1] 215   3

## [1] "q0002_hhid" "hea2"       "dep2"

## [1] 129   3

## [1] "q0002_hhid" "hea3"       "dep3"

Para poder combinarlas verticalmente, los nombres de las columnas que queremos combinar tienen que ser iguales. Para ello usamos otra función de la librería dplyr: rename

acr1 <- ac1 %>% # partimos de ac1 y entonces
  rename(hea = hea1) # renombramos hea1 como hea
acr2 <- ac2 %>% # partimos de ac2 y entonces
  rename(hea = hea2, dep = dep2) # renombramos hea2 como hea y dep2 como dea
acr3 <- ac3 %>% 
  rename(hea = hea3, dep = dep3) # mismo proceso para ac3

#names(ac1); names(ac2); names(ac3)

Finalmente combinamos las filas con otra función de la librería dplyr: bind_rows

#?bind_rows
acv <- bind_rows(acr1, acr2, acr3, .id = "wave")

Otras opciones para combinar bases de datos

Ejercicio:

Ver qué variables tiene acv; para qué se añadió .id = "wave"?; ejecutar rbind(ac1,ac2,ac3)

Nota: La función rbind de R hace esencialmente lo mismo, pero necesita que las 3 bases de datos tengan exactamente las mismas columnas. En caso en que esto no ocurre, como el presente, da un error. Además, tanto en el caso de rbind como de bind_rows, conviene que las columnas por las que se combina (las de igual nombre) tengan la misma clase, pues en caso contrario pueden ocurrir errores o comportamiento extraños.

Nota: Las funciones cbind y bind_cols de R y dplyr respectivamente combinan por columnas. Se diferencian de un merge en que no hay una columna “común” por la que fusionar, sino que se añaden las columnas de los distintos objetos, tal como están ordenadas en cada uno de ellos. Además, todas las columnas tienen que tener el mismo número de elementos.

Fusionamos horizontalmente las 3 bbdd originales de los ensayos clínicos. En este caso no las podemos combinar porque tienen distinto número de filas:

ach <- ac1 %>%
  full_join(ac2, by = c("q0002_hhid")) %>%
  full_join(ac3, by = c("q0002_hhid"))
ach %>% 
  head()

## # A tibble: 6 × 11
##   q0002_hhid sex   age   mar_1 edu1          hea1 grups         hea2 dep2   hea3 dep3 
##        <dbl> <chr> <chr> <chr> <chr>        <dbl> <chr>        <dbl> <chr> <dbl> <chr>
## 1       2306 fem   50-64 No    Primary       45.9 Intervention  NA   <NA>   NA   <NA> 
## 2       3363 fem   65-79 No    Less primary  23.3 Intervention  24.6 No     NA   <NA> 
## 3       3228 fem   50-64 Yes   Less primary  49.1 Intervention  35.2 No     NA   <NA> 
## 4       2371 fem   18-34 No    Secondary     57.5 Intervention  NA   <NA>   NA   <NA> 
## 5       2452 fem   35-49 No    Less primary  41.3 Intervention  44.3 No     31.6 Yes  
## 6       2750 fem   80+   No    Primary       19.5 Intervention  NA   <NA>   NA   <NA>

# cbind(ac1, ac2, ac3) %>% 
#    head()

Rellenando datos

Cuando se combinan las filas mediante bind_cols, las columnas que sólo están en una base de datos se llenan con missings. Por ejemplo, es el caso de las variables sociodemográficas, que en el dataframe acv sólo tienen datos en las filas correspondientes a la ola 1, que son las filas obtenidas de ac1.

acv %>%
  arrange(q0002_hhid) %>% # ordenamos por id
  head()

## # A tibble: 6 × 9
##   wave  q0002_hhid sex   age   mar_1 edu1        hea grups   dep  
##   <chr>      <dbl> <chr> <chr> <chr> <chr>     <dbl> <chr>   <chr>
## 1 1             48 fem   50-64 No    Secondary  47.9 Control <NA> 
## 2 2             48 <NA>  <NA>  <NA>  <NA>       46.3 <NA>    No   
## 3 3             48 <NA>  <NA>  <NA>  <NA>       47.4 <NA>    No   
## 4 1             61 fem   35-49 No    Tertiary   63.2 Control <NA> 
## 5 2             61 <NA>  <NA>  <NA>  <NA>       53.2 <NA>    No   
## 6 3             61 <NA>  <NA>  <NA>  <NA>       68.0 <NA>    No

Acudimos ahora a la función fill de la librería tidyr:

library(tidyr)
acv %>%
  group_by(q0002_hhid) %>% # agrupamos por id y entonces
  arrange(wave) %>% # ordenamos por ola para tener la fila con datos antes que las otras, y entonces
  fill(sex, age, mar_1, edu1, grups) %>% # para cada id rellenamos las filas vacías de las variables especificadas con los valores que no son missings, y entonces
  ungroup() %>% # desagrupamos y entonces
  arrange(q0002_hhid) %>% # ordenamos por id para ver las mismas filas que arriba, y entonces
  head() # nos quedamos con las primeras filas

## # A tibble: 6 × 9
##   wave  q0002_hhid sex   age   mar_1 edu1        hea grups   dep  
##   <chr>      <dbl> <chr> <chr> <chr> <chr>     <dbl> <chr>   <chr>
## 1 1             48 fem   50-64 No    Secondary  47.9 Control <NA> 
## 2 2             48 fem   50-64 No    Secondary  46.3 Control No   
## 3 3             48 fem   50-64 No    Secondary  47.4 Control No   
## 4 1             61 fem   35-49 No    Tertiary   63.2 Control <NA> 
## 5 2             61 fem   35-49 No    Tertiary   53.2 Control No   
## 6 3             61 fem   35-49 No    Tertiary   68.0 Control No

Pivotar (vertical/largo –> horizontal/ancho)

Pivotar es el proceso de transformar una base de datos vertical/longitudinal (por ejemplo el resultado de combinar filas, como es el caso de acv) a una horizontal (por ejemplo el producto de un una fusión, como es el caso de ach) o a la inversa.

Ejercicio:

• Modificar acv de tal manera que sus variables sociodemográficas estén completas tal como se indicó en la diapositiva anterior.
• Qué variables tiene acv?
• Qué variables de avc no tienen necesariamente los mismos valores en las distintas olas?

Pivotaje usando la función pivot_wider de la librería tidyr:

#?pivot_wider
acv %>%
  pivot_wider(names_from = c("wave"), values_from = c("hea", "dep")) %>%
  # select(-where(~all(is.na(.)))) %>% # quitamos aquellas columnas que cumplen que todos sus valores son missings
  head()

## # A tibble: 6 × 12
## # Groups:   q0002_hhid [6]
##   q0002_hhid sex   age   mar_1 edu1         grups        hea_1 hea_2 hea_3 dep_1 dep_2 dep_3
##        <dbl> <chr> <chr> <chr> <chr>        <chr>        <dbl> <dbl> <dbl> <chr> <chr> <chr>
## 1       2306 fem   50-64 No    Primary      Intervention  45.9  NA    NA   <NA>  <NA>  <NA> 
## 2       3363 fem   65-79 No    Less primary Intervention  23.3  24.6  NA   <NA>  No    <NA> 
## 3       3228 fem   50-64 Yes   Less primary Intervention  49.1  35.2  NA   <NA>  No    <NA> 
## 4       2371 fem   18-34 No    Secondary    Intervention  57.5  NA    NA   <NA>  <NA>  <NA> 
## 5       2452 fem   35-49 No    Less primary Intervention  41.3  44.3  31.6 <NA>  No    Yes  
## 6       2750 fem   80+   No    Primary      Intervention  19.5  NA    NA   <NA>  <NA>  <NA>

Para el pivotaje inverso se utilizaría la función pivot_longer.

#?pivot_longer

T-test

Comparamos entre grupos mediante un t-test:

t.test(hea1 ~ sex, data = dataw1)

## 
##  Welch Two Sample t-test
## 
## data:  hea1 by sex
## t = -14.557, df = 4523.6, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group fem and group masc is not equal to 0
## 95 percent confidence interval:
##  -7.861186 -5.995119
## sample estimates:
##  mean in group fem mean in group masc 
##           50.82465           57.75280

t.test(data$hea1, data$hea2, paired = TRUE) #comparamos el estado de salud en la ola 1 con la ola 2

## 
##  Paired t-test
## 
## data:  data$hea1 and data$hea2
## t = 0.81332, df = 2384, p-value = 0.4161
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.2407952  0.5820965
## sample estimates:
## mean of the differences 
##               0.1706506

Chi-cuadrado y correlación

-Test de chi-cuadrado

chisq.test(dataw1$dep1, dataw1$sex)

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  dataw1$dep1 and dataw1$sex
## X-squared = 68.537, df = 1, p-value < 2.2e-16

-Cálculo de correlaciones

#cor(dataw1$score_sup, dataw1$income)
cor.test(dataw1$score_sup1, dataw1$income1)

## 
##  Pearson's product-moment correlation
## 
## data:  dataw1$score_sup1 and dataw1$income1
## t = 3.2915, df = 3518, p-value = 0.001007
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.02241180 0.08828366
## sample estimates:
##        cor 
## 0.05540802

Ejercicio:

Consultar la ayuda de cor. Calcular la correlación de Spearman. Cómo tendríamos que hacer si hay valores missing?

Regresión lineal:

#?lm
fit <- lm(hea1 ~ age + sex + income1 + dep1, data = dataw1)
summary(fit)

## 
## Call:
## lm(formula = hea1 ~ age + sex + income1 + dep1, data = dataw1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -53.560  -8.844   0.710   9.444  40.641 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  64.99142    0.78354  82.946  < 2e-16 ***
## age35-49     -5.41946    0.95045  -5.702 1.28e-08 ***
## age50-64    -11.26679    0.79368 -14.196  < 2e-16 ***
## age65-79    -19.19973    0.80852 -23.747  < 2e-16 ***
## age80+      -30.15751    1.02706 -29.363  < 2e-16 ***
## sexmasc       5.07076    0.44117  11.494  < 2e-16 ***
## income1       0.16659    0.02742   6.075 1.37e-09 ***
## dep1Yes     -14.52514    0.67309 -21.580  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 13.09 on 3667 degrees of freedom
##   (1078 observations deleted due to missingness)
## Multiple R-squared:  0.3837, Adjusted R-squared:  0.3825 
## F-statistic: 326.1 on 7 and 3667 DF,  p-value: < 2.2e-16

Regresión logística

#?glm
fit <- glm(dep1 ~ age + sex + income1*hea1 + mar_1*score_lon1, data = dataw1, family = binomial)
summary(fit)

## 
## Call:
## glm(formula = dep1 ~ age + sex + income1 * hea1 + mar_1 * score_lon1, 
##     family = binomial, data = dataw1)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.0231  -0.4485  -0.2799  -0.1678   3.2085  
## 
## Coefficients:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)          0.9027594  0.4760209   1.896 0.057898 .  
## age35-49             0.1435643  0.3101981   0.463 0.643497    
## age50-64            -0.2696730  0.2772296  -0.973 0.330681    
## age65-79            -1.0277636  0.2922478  -3.517 0.000437 ***
## age80+              -2.0710935  0.3598878  -5.755 8.67e-09 ***
## sexmasc             -0.2519979  0.1314381  -1.917 0.055208 .  
## income1             -0.0154454  0.0276949  -0.558 0.577051    
## hea1                -0.0720689  0.0064107 -11.242  < 2e-16 ***
## mar_1Yes            -0.5831916  0.3266439  -1.785 0.074196 .  
## score_lon1           0.3268351  0.0406748   8.035 9.33e-16 ***
## income1:hea1        -0.0001597  0.0005689  -0.281 0.778968    
## mar_1Yes:score_lon1  0.1290971  0.0668269   1.932 0.053383 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 2585.4  on 3566  degrees of freedom
## Residual deviance: 1937.0  on 3555  degrees of freedom
##   (1186 observations deleted due to missingness)
## AIC: 1961
## 
## Number of Fisher Scoring iterations: 6

Intervalos de confianza

Con la función t.test podemos obtener el intervalo de confianza de una colección de valores

t.test(dataw1$hea1)$conf.int

## [1] 53.48650 54.44456
## attr(,"conf.level")
## [1] 0.95

Intervalos de confianza de los coeficientes de un modelo:

#?confint
confint(fit)

## Waiting for profiling to be done...

##                            2.5 %        97.5 %
## (Intercept)         -0.038058968  1.8302067963
## age35-49            -0.453392497  0.7675523589
## age50-64            -0.794681232  0.2966808469
## age65-79            -1.585558738 -0.4353302342
## age80+              -2.774032729 -1.3587725617
## sexmasc             -0.511606534  0.0040427621
## income1             -0.070008567  0.0387003866
## hea1                -0.084808502 -0.0596656506
## mar_1Yes            -1.223548026  0.0578792815
## score_lon1           0.247665974  0.4072594365
## income1:hea1        -0.001288554  0.0009438818
## mar_1Yes:score_lon1 -0.001686246  0.2604890641

Gràficas de modelos lineales simples usando ggplot: una única variable predictora

library(ggplot2)
fit <- lm(hea1 ~ income1, data = dataw1)
dataw1 %>%
  filter(!is.na(hea1), !is.na(income1)) %>% # nos quedamos con las observaciones que no tienen missings en el modelo
  # filter(!if_any(c(hea1, income1), is.na)) %>% # lo mismo de forma resumida
  cbind(predict(fit, interval = "confidence")) %>% # añadimos la predicción del modelo y los intervalos de confianza como nuevas columnas
  ggplot(aes(x = income1)) + # creamos el cuadro y añadimos el eje x
  geom_point(aes(y = hea1)) + # añadimos los punos d la base de datos
  geom_line(aes(y = fit)) + # añadimos la línea que el modelo induce
  geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = .15) + # añadimos los CI; también es posible como sigue
  # geom_line(aes(y = lwr), color = "red", linetype = "dashed") +
  # geom_line(aes(y = upr), color = "red", linetype = "dashed") +
  theme_bw() # cambiamos el tema por defecto

Gràficas de modelos lineales usando ggplot: una variable predictora continua y otra categórica

library(ggplot2)
fit <- lm(hea1 ~ income1 + edu1, data = dataw1)
dataw1 %>%
  # filter(!is.na(hea1), !is.na(income1), !is.na(edu1)) %>% # nos quedamos con las observaciones que no tienen missings en el modelo
  filter(!if_any(c(hea1, income1, edu1), is.na)) %>% # lo mismo de forma resumida
  cbind(predict(fit, interval = "confidence")) %>% # añadimos la predicción del modelo y los intervalos de confianza como nuevas columnas
  ggplot(aes(x = income1, color = edu1)) + # creamos el cuadro y añadimos el eje x
  geom_point(aes(y = hea1)) + # añadimos los punos d la base de datos
  geom_line(aes(y = fit)) + # añadimos la línea que el modelo induce
  geom_ribbon(aes(ymin = lwr, ymax = upr, color = NULL, fill = edu1), alpha = .15) # añadimos los CI

Gràficas de modelos lineales generalizados simples usando ggplot: una única variable predictora

library(ggplot2)
fit <- glm(dep1 ~ income1, data = dataw1, family = binomial)
dataw1 %>%
  filter(!if_any(c(dep1, income1), is.na)) %>%  # nos quedamos con las observaciones que no tienen missings en el modelo
  mutate(fit = predict(fit, type = "response")) %>% # añadimos la predicción del modelo como nuevas columnas
  ggplot(aes(x = income1)) + # creamos el cuadro y añadimos el eje x
  geom_line(aes(y = fit)) + # añadimos la línea que el modelo induce
  labs(y = "Predicted probability of dep1 == Yes") + # especificamos la etiqueta del eje y
  theme_classic() #cambiamos el tema por defecto

Gràficas de modelos lineales generalizados usando ggplot: una variable predictora continua y otra categórica

library(ggplot2)
fit <- glm(dep1 ~ income1 + edu1, data = dataw1, family = binomial)
dataw1 %>%
  filter(!if_any(c(dep1, income1, edu1), is.na)) %>%  # nos quedamos con las observaciones que no tienen missings en el modelo
  mutate(fit = predict(fit, type = "response")) %>% # añadimos la predicción del modelo como nuevas columnas
  ggplot(aes(x = income1, color = edu1)) + # creamos el cuadro y añadimos el eje x
  geom_line(aes(y = fit)) + # añadimos la línea que el modelo induce
  labs(y = "Predicted probability of dep1 == Yes") + # especificamos la etiqueta del eje y
  theme_minimal() # cambiamos el tema por defecto

Otras formas de obtener resultados estadísticos en R

De la misma manera que tidyverse hace referencia a una colección de librerías que simplifican la manipulación de datframes, existe otra colección de librerías, easystats (https://easystats.github.io/easystats/), cuyo objetivo es mostrar los resultados de tests y modelos estadísticos con mejor formato.

Ejercicio 8:

Instalar y cargar las librerías performance y parameters y probar las funciones model_performance y model_parameters aplicadas a fit.

La librería modelsummary dispone de 2 funciones que permiten visualizar multiples resultados en tablas bien organizadas:

• modelsummary (https://vincentarelbundock.github.io/modelsummary/articles/modelsummary.html) y
• datasummary (https://vincentarelbundock.github.io/modelsummary/articles/datasummary.html).

Otras librerías que permiten visualizar descriptivos y resultados resumidos en tablas son:

• summarytools (https://cran.r-project.org/web/packages/summarytools/vignettes/introduction.html) y
• compareGroups (https://cran.r-project.org/web/packages/compareGroups/vignettes/compareGroups_vignette.html).

Existen distintas librerías que facilitan la representación gráfica de modelos, predicciones, …:

• jtools (https://cran.r-project.org/web/packages/jtools/vignettes/effect_plot.html)
• effects (https://cran.r-project.org/web/packages/effects/vignettes/predictor-effects-gallery.pdf)
• ggeffects(https://strengejacke.github.io/ggeffects/)

Algunos recursos

Tutoriales introductorios a R:

• http://rosuda.org/mitarbeiter/pilhoefer/rkurs2.pdf
• http://b2slab.upc.edu/software-and-tutorials/r-nutshell/
• https://www.cyclismo.org/tutorial/R/index.html
• https://www.uv.es/vcoll/preliminares.html
• http://people.math.aau.dk/~sorenh/misc/rdocs/Rintro-notes.pdf
• https://cran.r-project.org/doc/contrib/Verzani-SimpleR.pdf
• en la pestaña Resources de la web de RStudio

Blogs:

• https://statisticsglobe.com/r-programming-language/
• enfocados a la estadística, pero no sólo (de hecho, diría que los siguientes son los mejores):
  • http://www.sthda.com/english/
  • http://www.flutterbys.com.au/stats/course.html
  • http://www.r-tutor.com/ (incluye contenido de estadística bayesiana y machine learning, por otra parte un poco anticuado)
• enfocado a la investigación en psicología (y vinculado a la librería psych):
  • http://personality-project.org/r/
• https://www.statmethods.net/index.html

Multitud de libros:

• https://bookdown.org/
• https://www.bigbookofr.com/ (contenidos parecidos a los del link anterior, en menor cantidad, aunque clasificados por materias)

Cursos online (MOOC’s) en plataformas como Coursera, edX, etc.

Más recursos

Acerca de librerías o materias específicas:

• sobre librerías en general:
  • https://rdrr.io/
  • https://cran.r-project.org/web/views/
• gráficas en general:
  • https://www.r-graph-gallery.com/
  • https://r-charts.com
• formas de mostrar descriptivos:
  • https://dabblingwithdata.wordpress.com/2018/01/02/my-favourite-r-package-for-summarising-data/
• meta-análisis:
  • https://bookdown.org/MathiasHarrer/Doing_Meta_Analysis_in_R/
• análisis factorial:

Recursos de RStudio

• El foro https://community.rstudio.com/ principalmente centrado (aunque no exclusivamente) en cuestiones relacionadas con RStudio y con las librerías de la familia tidyverse
• La página web https://support.rstudio.com/hc/en-us contiene material acerca de RStudio
• El blog https://blog.rstudio.com publica contenido diverso relacionado con los productos de la compañía RStudio
• El blog https://rviews.rstudio.com está más enfocado al usuario de R, aportando información y ejemplos del uso de distintas librerías, así como acerca de conferencias de usuarios de R
• El blog https://www.tidyverse.org/blog/ publica habitualmente entradas acerca de las nuevas librerías de la familia tidyverse o de actualizaciones relevantes